Tight— binding description of the quasiparticle dispersion of graphite and few— layer 

graphene 



A. Griineis 1 * 2 *, C. Attaccalitc 3 ,L. Wirtz 4 , H. Shiozawa 5 , R. Saito 6 ,T. Pichler 1 , A. Rubio 3 
1 Faculty of Physics, Vienna University, Strudlhofgasse 4, 1090 Wien, Austria 
2 IFW Dresden, P.O. Box 270116, D-01171 Dresden, Germany 
3 Dept. Fisica de Materiales, Donostia International Physics Center, Spain 
European Theoretical Spectroscopy Facility (ETSF), Spain 
4 Institute for Electronics, Microelectronics, and Nanotechnology Dept. ISEN B.P. 60069 59652 Villeneuve d'Ascq Cedex, France 
5 Advanced Technology Institute, University of Surrey, Guildford, GU2 7XH, UK and 
6 Department of Physics, Tohoku University, Aoba, Sendai, 980-8578, Japan 
(Dated: August 13, 2008) 

A universal set of third-nearest neighbour tight-binding (TB) parameters is presented for calcu- 
lation of the quasiparticle (QP) dispersion of N stacked sp 2 graphene layers (N = 1 . . . oo) with AB 
stacking sequence. The QP bands are strongly renormalized by electron-electron interactions which 
results in a 20% increase of the nearest neighbour in-plane and out-of-plane TB parameters when 
compared to band structure from density functional theory. With the new set of TB parameters 
we determine the Fermi surface and evaluate exciton energies, charge carrier plasmon frequencies 
and the conductivities which are relevant for recent angle-resolved photoemission, optical, electron 
energy loss and transport measurements. A comparision of these quantitities to experiments yields 
an excellent agreement. Furthermore we discuss the transition from few layer graphene to graphite 
and a semimetal to metal transition in a TB framework. 



1. INTRODUCTION 

Recently mono- and few-layer graphene (FLG) in an 
AB (or Bernal) stacking is made with high crystallinity 
by the following three methods; epitaxial growth on 
SiC PQ, chemical vapour deposition on Ni(lll) [2] and 
by mechanical cleavage on SiC>2 [3J ■ Graphene is a novel, 
two-dimensional (2D) and meta stable material which 
has sparked interest from both basic science and appli- 
cation point of view [4]. A monolayer of graphene al- 
lows one to treat basic questions of quantum mechanics 
such as Dirac Fermions or the Klein paradox ^ in a 
simple condensed-matter experiment. The existence of 
a tunable gap in a graphene bilayer was shown by angle- 
resolved photoemission (ARPES) jS], which offers a pos- 
sibility of using these materials as transistors in future 
nanoelectronic devices that can be lithographically pat- 
terned [7]. Furthermore a graphene layer that is grown 
epitaxially on a Ni(lll) surface is a perfect spin filter 
device [5] that might find applications in organic spin- 
tronics. 

It was shown recently by ARPES that the electronic 
structure of graphene j9] and its 3D parent material, 
graphite [10] [IT] [12], is strongly renormalized by cor- 
relation effects. To date the best agreement between 
ARPES and ab-initio calculations is obtained for GW 
(Greens function G of the Coulomb interaction W) calcu- 
lations of the QP dispersion. The band structure in the 
local density approximation (LDA) (bare band disper- 
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sion) is not in good agreement with the ARPES spectra 
because it does not include long-range correlation effects. 
The self energy correction of the Coulomb interaction to 
the bare energy band structure are crucial for determin- 
ing the transport and optical properties (excitons) and 
related condensed-matter phenomena. For graphite, a 
semi metal with a tiny Fermi surface, the number of 
free electrons to screen the Coulomb interaction is low 
(~ 10 19 carriers cm -3 ) and thus the electron-electron 
correlation is a major contribution to the self-energy 
correction [10]. Theoretically the bare energy band dis- 
persion is calculated by the local density approximation 
(LDA) and the interacting QP dispersion is obtained by 
the GW approximation. The GW calculations are com- 
putationally expensive and thus only selected k points 
have been calculated [TU]. Therefore a tight-binding 
(TB) Hamiltonian with a transferable set of TB param- 
eters that reproduces the QP dispersion in sp 2 stacked 
graphene sheets is needed for analysis of ARPES, optical 
spectroscopies and transport properties for pristine and 
doped graphite and FLGs. So far there are already sev- 
eral sets of TB parameters published for graphene, FLG 
and graphite. For graphene a third nearest neighbour fit 
to LDA has been performed [13] • Recently, however it 
has been shoen by ARPES that the LDA underestimates 
the slope of the bands and also the trigonal warping ef- 
fect [TU]. For bilayer graphene the parameters of the 
so-called Slonzcewski-Weiss-McClure (SWMC) Hamil- 
tonian have been fitted to reproduce double resonance 
Raman data [14]. A direct observation of the quasipar- 
ticle (QP) band structure is possible by ARPES. A set 
of TB parameters has been fitted to the experimental 
ARPES data of graphene grown on SiC [IB] ■ As a result 
they obtained a surprisingly large absolute value of the 
nearest neighbour ir hopping parameter of 5.13 eV |15j . 
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This is in stark contrast to the fit to the LDA calculation 
which gives only about half of this value [T3]. Consid- 
ering the wide range of values reported for the hopping 
parameters, a reliable and universal set of TB parame- 
ters is needed that can be used to calculate the QP dis- 
persion of an arbitrary number of graphene layers. The 
band structure of FLGs has been calculated [16 using 
first-nearest neighbour in-plane coupling which provide 
the correct band structure close to K point. However, 
as we will show in detail in this paper, the inclusion of 
third-nearest neighbours (3NN) is essential in describ- 
ing the experimental band structure in the whole BZ. 
The fact that the inclusion of 3NN is essential is also 
proven by the dispersion of a localized state at the zig- 
zag edge of a graphene flake. Only inclusion of 3NN in- 
teraction can reproduce a weakly downwards dispersing 
state which is relevant to superconductivity in graphene 
nanoribbons [T?] . 

In this paper, we present a tight -binding (TB) formu- 
lation of the 7r bare energy band and QP dispersions of 
AB stacked FLG and graphite. We have previously com- 
pared both, GW and LDA calculations to ARPES exper- 
iments and proofed that LDA underestimates the slope of 
the bands and trigonal warping [T5]. Here we list the TB 
fit parameters of the QP dispersion (TB-GW) and the 
bare band dispersion (TB-LDA) and show that the in- 
plane and out-of-plane hopping parameter increase when 
going from LDA to GW. This new and improved TB- 
GW parameters is used for direct comparision to exper- 
iments are obtained from a fit to QP calculations in the 
GW approximation. This set of TB parameters works in 
the whole 2D (3D) BZ of FLG (graphite) and is in agree- 
ment to recent experiments. In addition we fit the pa- 
rameters of the popular SWMC Hamiltonian that is valid 
close to the KH axis of graphite. This paper is organized 
as follows: in section 2 we develop the 3NN TB formula- 
tion for graphite and FLGs and in section 3 the SWMC 
Hamiltonian is revised. In section 4 a new set of TB pa- 
rameters for the calculating the QP dispersion of stacked 
sp 2 carbon is given. In section 5 we compare the graphite 
bare energy band (LDA) to the QP (GW) dispersion. In 
section 6 we use the TB-GW Hamiltonian and calculate 
the doping dependent Fermi surface of graphite and es- 
timate effective masses and free charge carrier plasmon 
frequencies of pristine and doped graphite. In section 7 
we show the calculated QP dispersions of FLGs. In sec- 
tion 8 we discuss the present results and estimate the ex- 
citon binding energies, transport properties and the low 
energy plasmon frequencies. In section 9 the conclusions 
of this work are given. Finally, in the appendices, the 
analytical forms of the Hamiltonians for FLG are shown. 



2. THIRD-NEAREST NEIGHBOUR TIGHT 
BINDING FORMULATION 

Natural graphite occurs mainly with AB stacking or- 
der and has four atoms in the unit cell (two atoms for 




FIG. 1: (a) The graphite unit cell consists of four atoms de- 
noted by Ai,Bi,A2 and B2 (light blue). The red arrows de- 
note the interatomic tight-binding hopping matrix elements 
7o, 70, 70, 71, • • • , 75. The overlap matrix elements so — S3 
(npot shown) couple the same atoms as 70 — 70 • (b) The 3D 
Brillouin zone of graphite with the high symmetry points and 
the coordinate system used throughout this work. 



each graphene plane) as shown in Fig. [TJa). Each atom 
contributes one electron to the four 7r electronic energy 
bands in the 3D Brillouin zone (BZ) [see Fig. [IJb)]. 
FLG has N parallel graphene planes stacked in an AB 
fashion above one another; the unit cell of FLG is 2D 
and the number of ir bands in the 2D BZ equals 27V. 
For graphite and FLG the TB calculations are carried 
out with a new 3NN Hamiltonian and in addition with 
the well-known SWMC Hamiltonian [T5] that is valid 
in the vicinity to the Fermi level (Ep). The TB pa- 
rameters that enter these two Hamiltonians are 7 = 
(7o > 7o > 7o > s o, si, s 3 , 7i • • • > 75, A, E ) for the 3NN Hamil- 
tonian [shown in Fig. l|a)] and 7' = (70, . . . , 75, A', E' ) 
for the SWMC Hamiltonian. The hopping matrix ele- 
ments for the SWMC Hamiltonian are not shown here 
but they have a similar meaning with the difference that 
only one nearest neighbour in-plane coupling constant 
is considered (see e.g. Ref. [TS] for an explanation of 
SWMC parameters). The hopping matrix elements for 
the 3NN Hamiltonian are shown in Fig.Qa). The atoms 
in the 3D unit cell are labelled Ax,Bi for the first layer 
and ^2,-82 for the second layer. The A2 atom lies directly 
above the A\ atom in z direction (perpendicular to the 
layers). Within the xy plane the interaction is described 
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by 7q (e.g. AiBi and A 2 B 2 ) for the nearest neighbours, 
and 7q and 7q for second nearest and third nearest neigh- 
bours, respectively. A further parameter, 71 (A\B 2 ), is 
needed to couple the atoms directly above each other 
(in z direction). The hopping between adjacent layers of 
sites that do not lie directly above each other is described 
by 73 (B1B2) and 74 (A1B2 and BiA 2 ). The small cou- 
pling of atoms in the next-nearest layer is 72 (B\B\ and 
B 2 B 2 ) and 75 {AxAt and A 2 A 2 ). 

The calculation shown here is valid for both graphite 
and FLGs with small adjustment as indicated when 
needed. The lattice vectors for graphite in the xy plane 
are ai and a% and the out-of-plane lattice vector per- 
pendicular to the layers is a.3. 

ai = (^,f,0), aa = (^ I -* 0), a 3 = (0,0,2c ) 

The C-C distance ao = 1.42A and the distance of two 
graphene layers Co = 3.35A. For a FLG with N layers 
and hence 2N atoms only the 2D unit vectors ai and 
a.2- Similarly the electron wave vectors in graphite k = 
(k x , k y , k z ) have 3 components and in FLG k = (k x , k y ). 
A TB method (or linear combination of atomic orbitals, 
LCAO) is used to calculate the bare energy band and QP 
dispersion by two different sets of interatomic hopping 
matrix elements. The electronic cigcnfunction ^(r,k) is 
made up from a linear combination of atomic 2p z orbitals 
cj)(r) which form the tt electronic bands in the solid. The 
electron wave function for the band with index j is given 

by 

*F(k,r) = Yl ci(k)<Mk,r), (j = 1...2N). 

s=Ai,Bi,...,B N 

Here j = 1 . . . 2N is the electronic energy band index and 
s in the sum is taken over all atomic 2p z orbitals from 
atoms Ax,Bi,A 2 ■ ■ ■ Bis; . Note that for 3D graphite we 
have N = 2. The c 3 s (k) are wave function coefficients for 
the Bloch functions $ s (k,r). The Bloch wave functions 
are given by a sum over the atomic wave functions </> s for 
each orbital in the unit cell with index (£, m, m) multi- 
plied by a phase factor. The Bloch function in graphite 
for the atom with index s is given by 

1 U 

$ s (r,k) = — £ exp(- J Rf m "-k)0 s (r-Rf m "). (3) 

where U is the number of unit cells and <f) s denotes atomic 
wave functions of orbital s. In graphite the atomic orbital 
<t> s in the unit cell with index (£, m, m) is centered at 
Rf mn = &li + ma 2 + na 3 + r s and in FLGs at R l s m = 
£ai + ma 2 + r s . 

For the case of FLGs <& s (r,k) contains only the sum 
over the 2D in-plane Rf™. The 2N x 2N Hamiltonian 
matrix defined by H ss /(k) = ($ s (r,k)|J?(r)|$ s /(r,k)) 
and the overlap matrix is defined by S ss / (k) = 
($ s (r, k)|$ s /(r, k)). For calculation of H and S, up to 



third nearest neighbour interactions (in the xy plane) and 
both nearest and next-nearest neighbour planes (in z di- 
rection) are included as shown in Fig. [l|a). The energy 
dispersion relations are given by the eigenvalues E(k) 
and are calculated by solving 

#(k)c(k) = 5(k)E(k)c(k). (4) 

In the Appendix, we show the explicit form of H (k) for 
graphite FLGs with 1-3 layers. 

3. SWMC HAMILTONIAN 

The SWMC Hamiltonian has been extensively used 
in the literature [19]. It considers only first nearest 
neighbour hopping and is valid close to the KH axis of 
graphite. For small k measured from the KH axis (up 
to 0.15 A- 1 ) both, the 3NN and the SWMC Hamiltoni- 
ans yield identical results. The eight TB parameters for 
the SWMC Hamiltonian were previously fitted to various 
optical and transport experiments |19j . 

For the cross sections of the electron and hole pock- 
ets analytical solutions have been obtained and thus it 
has been used to calculate the electronic transport prop- 
erties of graphite. Thus, in order to provide a connec- 
tion to many transport experiments from the past, we 
also fitted the LDA and GW calculations to the SWMC 
Hamiltonian. 

The TB parameters are directly related to the energy 
band structure. E.g. 70 is proportional to the Fermi 
velocity in the k x k y plane and 471 gives the bandwidth 
in the k z direction. The bandwidth of a weakly disper- 
sive band in k z that crosses Ep approximately halfway 
in between K and H is equal to 272, which is responsi- 
ble for the semi metallic character of graphite. Its sign 
is of great importance for the location of the electron 
and hole pockets: a negative sign brings the electron 
pocket to K while a positive sign brings the electron 
pocket to H. There has been positive signs of 72 re- 
ported earlier [5D] but it has been found by M.S. Dressel- 
haus [TSJ |2T] that the electron (hole) pockets are located 
at K(H) which is in agreement with recent DFT cal- 
culations [ID], tight-binding calculation [22] and experi- 
ments [HO [53] • The magnitude of 72 determines the over- 
lap of electrons and holes and the volume of the Fermi 
surface. It thus also strongly affects the concentration 
of carriers and hence the conductivity and free charge 
carrier plasmon frequency. The effective masses for elec- 
trons and holes of the weakly dispersive energy band are 
denoted by m* e and m* zhl respectively. Their huge value 
also results from the small value of 72 and causes the 
low electrical conductivity and low plasmon frequency in 
the direction perpendicular to the graphene layers since 
m ze(V m ze) enters the denominator in the expression for 
the Drude conductivity (free carrier plasmon frequency) . 
73 determines the strength of the trigonal warping effect 
(73 = gives isotropic equi-energy contours) and 74 the 
asymmetry of the effective masses in valence band (VB) 
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Method 


70 1 


7o 2 


7o 3 


so 


Sl 


S2 


7i 


72 


73 


74 


75 


E 


A 


3NN TB-GW 

3NN TB-LDA 


-3.4416 
-3.0121 


-0.7544 
-0.6346 


-0.4246 
-0.3628 


0.2671 
0.2499 


0.0494 
0.0390 


0.0345 
0.0322 


0.3513 
0.3077 


-0.0105 
-0.0077 


0.2973 
0.2583 


0.1954 
0.1735 


0.0187 
0.0147 


-2.2624 
-1.9037 


0.0540° 
0.0214 



"We adjusted the impurity doping level in order to reproduce the 
experimental value of A. 

TABLE I: 3NN tight-binding parameters for few-layer graphene and graphite. The parameters of fits to LDA and GW 
calculations are shown. The 3NN Hamiltonian is valid in the whole two(three) dimensional BZ of graphite (graphene layers). 
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■y'o 


7i 


72 
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74 


75 


K 


A' 


TB-GW" 


3.053 


0.403 


-0.025 


0.274 


0.143 


0.030 


-0.025 


-coos 6 


TB-LDA" 


2.553 


0.343 


-0.018 


0.180 


0.173 


0.018 


-0.022 


-0.018 


EXP C 


3.16 


0.39 


-0.02 


0.315 


0.044 


0.038 


-0.024 


-0.008 


LDA d 


2.598 


0.364 


-0.014 


0.319 


0.177 


0.036 


-0.026 


-0.013 


EXP e 


2.9 


0.3 




0.1 


0.12 








kkr/ 


2.92 


0.27 


-0.022 


0.15 


0.10 


0.0063 


0.0079 


-0.027 



"This work 

'We adjusted the impurity doping level to reproduce the experi- 
mental value of A' 
c Fit to Experiment, M.S. Dresselhaus et al 1191 
d Fit to LDA, J.C. Charlier et al [22] 

c Fit to double resonance Raman spectra. L.M. Malard et al. I14| 
^Fit to Korringa-Kohn-Rostocker first principles calculation 
Tatar and Rabi [18] 

TABLE II: The SWMC tight-binding parameters for the 
bare band dispersion (LDA) and the quasiparticle dispersion 
(GW). All values are in eV. This parameters are for the 
SWMC Hamiltonian [TO] OH E3 US] which is valid close to 
the KH axis. 



and conduction band (CB). The other parameter from 
next nearest neighbour coupling, 75 has less impact on 
the electronic structure: both the VB and CB at K are 
shifted with respect to the Fermi level by A + 75 causing 
a small asymmetry [IS]. Here A is the difference in the 
on-site potentials at sites A\(B2) and ^(.Bi). A is the 
value of the gap at the H point [24 . This crystal field 
effect for nonzero A occurs in AB stacked graphite and 
FLGs but it does not occur in AA stacked graphite and 
the graphene monolayer. The small on-site energy dif- 
ference A appears in the diagonal elements of i?(k). It 
causes an opening of a gap at the H point which results 
in a breakdown of Dirac Fermions in graphite and FLGs 
with N > 1. Finally, Eq is set in such a way that the 
electron and hole like Fermi surfaces of graphite yield an 
equal number of free carriers. Eq is measured from the 
bottom of the CB to the Fermi level. 



4. NUMERICAL FITTING PROCEDURE 

The ab-initio calculations of the electronic dispersion 
are performed on two levels: bare band dispersion cal- 
culation by LDA and QP dispersion calculations within 
the GW approximation. We calculate the Kohn-Sham 
band-structure within the LDA to density-functional the- 
ory (DFT) [26]. Wave- functions are expanded in plane 
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FIG. 2: The TB fit along k z direction at k z = k y = (KH 
axis) for (a) LDA and (b) GW. □ denotes LDA calculations 
and A denotes GW calculations taken from Ref. [10] that 
were used for the fitting. 



waves with an energy cutoff at 25 Ha. Core electrons are 
accounted for by Trouiller-Martins pseudopotentials. 

We then employ the GqWq approximation using a 
plasmon-pole approximation for the screening |27l |28j [29] 
to calculate the self-energy corrections to the LDA dis- 
persion. For the calculation of the dielectric function 
e(oj,q) we use a 15x15x5 Monkhorst-Pack k sampling of 
the first BZ, and conduction band states with energies 
up to 100 eV above the valence band (80 bands), cal- 
culations were performed using the code YAMBO |30j . 
The details for the first principles calculations are given 
elsewhere [5Tj . 

For the fitting of the TB parameters to the ab-initio 
(LDA and GW) calculations we used energies of the four 
7T bands of graphite at ~ 100 k points. The points were 
distributed inside the whole 3D BZ of graphite. The 
fitting was performed with the 3NN Hamiltonian. In 
addition we chose a smaller subset of points inside a 
volume of 0.15A- 1 x 0.15A- 1 x 0.47A- 1 and fitted the 
SWMC parameters, which is frequently used in the liter- 
ature [isniicre]). 

The set of TB parameters were fitted by employing 
a steepest-descent algorithm that minimizes the sum of 
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FIG. 3: Upper panel: The TB fits along k y direction at k z = 
(K point) for (a) LDA and (b) GW. Lower panel: the k y dis- 
persion at k z = 0.47A -1 (H point) for (c) LDA and (d) GW. 
□ denotes LDA calculations and A denotes GW calculations 
that were used for the fitting (taken from Ref. [10]). 



squared differences between the TB and the ab-initio cal- 
culations. This involves solving Eq. [4] with different sets 
of TB parameters so as to approach a minimum deviation 
from the ab-initio calculations. Points close to Ep were 
given additional weight so that the band crossing Ep was 
described with a deviation less than 1 meV. This is im- 
portant for an accurate description of the Fermi surface. 
In Table|T] we list the parameters for the 3NN Hamilto- 
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FIG. 4: The TB fits of equi-energy contours at fc z =0 for (a) 
TB-LDA and (b) TB-GW. 



nian that can be used to calculate TB-GW bands in the 
whole 3D BZ of graphite. 

The parameters that were fit with the SWMC Hamil- 
tonian are summarized in Table[TTJ These TB parame- 
ters reproduce the bare energy band (fit to LDA) and 
the QP (fit to GW) calculated dispersions. Hereafter 
these fits are referred to as TB-LDA and TB-GW, re- 
spectively. We also list the values from other groups 
that were fit to experiments [TH] and to LDA [22] and 
another first-principles calculation [18]. It can be seen 
that the TB-GW 7 parameters for the nearest neighbour 
coupling increase by about 20% when compared to TB- 
LDA. TB-GW is also closer to the experimental TB pa- 
rameters than the TB-LDA parameters. This indicates 
that electronic correlation effects play a crucial role in 
graphite and FLGs for interpreting and understanding 
experiments that probe the electronic energy band struc- 
ture. 



5. COMPARISION OF THE BARE ENERGY 
BAND TO THE QUASIPARTICLE DISPERSION 
OF GRAPHITE 



We now compare the calculated TB-LDA to TB-GW 7 
and we also show the result of the first-principles calcu- 
lations that were used for fitting in order to illustrate the 
quality of the fit. In Fig. [2j the full k z dispersion from K 
to H for (a) TB-LDA is compared to (b) TB-GW calcu- 
lations. It is clear that the bandwidth in k z increases by 
about 20% or 200 meV when going from TB-LDA to TB- 
GW, i.e. when long-range electron-electron interaction 
is taken into account. Such an increase in bandwidth is 
reflected by the TB parameter 71 (in the SWMC model 
471 is the total bandwidth in the out-of-plane direction) . 
It can be seen that the conduction bandwidth increases 
even more than the valence bandwidth. The VB disper- 
sion was measured directly by ARPES and gave a result 
in good agreement to the TB-GW [10] . 
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FIG. 5: The in-plane quasiparticle dispersion for (a) k z — and (b) fc z =0.47 A 1 calculated by TB-GW. The symbols 7Ti — 1:4 
denote the four n bands of graphite. 



The dispersion parallel to the layers is investigated in 
Fig. [3] where we show the k y dispersion for k z = ( K) 
and k z = 0.47A" 1 (H). Here the TB-GW bands also be- 
come steeper by about 20% when compared to TB-LDA. 
This affects 70 which is the in-plane nearest-neighbour 
coupling (see Fig. [T|(a). It determines the in-plane vf 
and in-plane bandwidth which is proportional to 7q — 7q 
(or proportional to 70 in the SWMC Hamiltonian) . In 
Fig.|4]the trigonal warping effect is illustrated by plotting 
an equi-energy contour for k z — with (a) TB-LDA and 
(b) TB-GW. The trigonal warping effect is determined 
by 73 which is larger in the TB-GW fit compared to the 
TB-LDA. 



6. THE THREE DIMENSIONAL 
QUASIPARTICLE DISPERSION AND DOPING 
DEPENDENT FERMI SURFACE OF GRAPHITE 



In the previous section we have shown that the ir band- 
width increases by 20% when going from TB-LDA to TB- 
GW . This affects especially the optical properties such 
as the 7r — > 7r* transition that plays an important role 
in optical absorption and resonance Raman and thus one 
has to use TB-GW for proper description of the elec- 
tronic structure of graphite including electron-electron 
correlation effects. In Fig. [5] we show the complete in- 
plane QP band structure for (a) k z — (K point) and 
(b) k z = 0.46 A" 1 (H point) calculated by the 3NN 
TB Hamiltonian. In the YKM plane two valence bands 
(•7Ti and 7T2) and two conduction bands (^3 and K4) can 
be seen and in the whole HAL plane the two valence 
(conduction) bands are degenerate. The ab-initio GW 
values are compared to the 3NN TB-GW calculation in 
Table III It can be seen that the fit reproduces the ab- 
initio calculations with an accuracy of ~ 10 meV in the 



Point 


Method 


7T1 


7T2 


7T3 


7T 4 


r 


GW 
TB-GW 


-9.458 
-9.457 


-7.257 
-7.258 


12.176 
12.184 


12.541 
12.540 


M 


GW 
TB-GW 


-3.232 
-3.216 


-2.441 
-2.457 


1.655 
1.656 


2.491 
2.495 


K 


GW 
TB-GW 


-0.736 
-0.728 


-0.025 
-0.024 


-0.025 
-0.024 


0.917 
0.909 


H 


GW 
TB-GW 


0.020 
0.020 


0.020 
0.020 


0.025 
0.025 


0.025 
0.025 



TABLE III: Energy values of the GW calculation and the 
3NN TB ht (TB-GW) at high symmetry points in the 3D 
BZ of graphite (all values in units eV). The symbols 7Ti-7T4 
denote the four n bands of graphite. The TB-GW have been 
calculated with the 3NN TB-GW parameters from Table I. 



BZ center and an accuracy of 1 meV along the KH axis, 
close to Ep. 

Close to Ep the 3NN Hamiltonian is identical to the 
SWMC Hamiltonian. Thus, for evaluation of the Fermi 
surface and the doping dependence on Ep in the dilute 
limit, we use the SWMC Hamiltonian. The weakly dis- 
persing band that crosses Ep is responsible for the Fermi 
surface and the electron and hole pockets. This energy 
band is illustrated in Fig{6]ja) where we show that the 
TB-GW fit has an accuracy of 1 meV. The minority 
pocket that is a result of the steeply dispersive energy 
band close to H is shown in Fig{6jb). 

The QP dispersion close to Ep is shown in FigjT] 
for (a,b) K, (c,d) H point. The dispersion around 
the K point is particular complicated: there arc four 
touching points between valence and conduction bands. 
Three touching points between the valence and conduc- 
tion bands exist in the close vicinity to K at angles of 
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FIG. 6: (a) The weakly dispersive band that is responsible 
for the formation of electron and hole pockets. A denotes 
GW calculations (with the impurity doping level adjusted to 
reproduce the experimental gap at H point) and the line the 
TB-GW fits, (b) shows a magnification of the minority hole 
pocket that is caused by the steeply dispersive band close to 
H. 



0°, 120° and 240° away from k x (i.e. the KM direction). 
The fourth touching point is exactly at K point. The 
touching points arise from the semi metallic character 
of graphite: there are two parabolas (VB and CB) that 
overlap by about 20 meV. For example the bottom of the 
CB is denoted by the point 2 in FigQa) and (b). At H 
point shown in Figj7]Jc) and (d) the energy band structure 
becomes simpler: there are only two non-degenerate en- 
ergy bands and their dispersion is rather isotropic around 
H (the trigonal warping effect in the k x k y plane is a min- 
imum in the AHL plane and a maximum in the TKM 
plane). It is clear that the energy bands do not touch 
each other and the dispersion is not linear but parabolic 
with a very large curvature (and hence a very small ab- 
solute value of the effective mass; see Fig. [8| at H. Ep 
lies 20 meV below the top of the VB and the energy gap 
A is equal to 5 meV. It is interesting to note that the a 
larger value of A would bring the top of the VB above 



Ep. The horizontal cuts through the dispersions in FigjT] 
at E — Ep (blue area) give cross sectional areas of an 
electron-like Fermi surface at k z = (K) and a hole-like 
Fermi surface at k z = 0.47 A -1 (H point), consistent 
with a semi-metallic behaviour. 

The effective in-plane massses for electrons (to*) and 
holes (to,*) arc evaluated along the parabola indicated in 
Fig|7]by the points 1,2,3. This parabola lies in the plane 
spanned by KM and HL and can thus be considered an 
upper limit for the effective mass since the dispersion is 
flat in this direction as can be seen in FigfTja). For to: 
the center of the parabola is chosen to be the bottom of 
the CB, i.e. point 2 in Figj7^a). Similarly for m* h the 
center of the parabola is chosen to be the top of the VB, 
i.e. point 2 in Fig{7jb). Due to the larger curvature of the 
hole bands, the absolute value of to* is larger than m* h . 
The k z dependence of to* and m* h is shown in Fig 8^a) 
and (b), respectively. For the effective mass in the z 
direction we fit a parabola for the weakly dispersing band 
in direction perpendicular to the layers and get m* e — 



16mn and 



'zh 



= — 16toq for the effective electron and 



hole masses perpendicular to the layers, respectively. The 
Figgc) and Fig^d) shows the TB-GW along k z for 
the electron and hole pocket, respectively. The parabolic 
fits that were used to determine the effective masses are 
shown along with the calculation. 

The k z dependence of to* and m* h are shown in Fig. [9] 
For this purpose we evaluated the heavy electron mass of 
the parabolic sub bands as shown in Fig. [7] and Fig. [8}. 
It is clear that to* has a weak k z dependence and m* h 
strongly depends on the value of k z . This is obvious since 
exactly at H point, the value of m* h has a minimum. For 
a finite value of the gap A , the value of m* h also remains 
finite. 

We now discuss the whole 3D Fermi surface. The vol- 
ume inside the surface determines the low-energy free 
carrier plasmon frequencies and the electrical conductiv- 
ity. The trigonal warping has little effect on the vol- 
ume inside the electron and hole pocket. When we 
set 73 = 0, then the Fermi surface is isotropic around 
KH axis. The simplification of 73 = results in lit- 
tle change of the volume. For 73 7^ 0, there are touch- 
ing points of the electron-like and hole-like Fermi sur- 
faces [IS]. The touching points (or legs) are important 
for understanding the period for de-Haas-van Alphcn 
and the large diamagnctism in graphite |25j . However, 
for the calculation of the number of carriers, they are 
not crucial and thus the Fermi surface calculated with 
73 = can be used for the evaluation of the electron 
density, n e , and the hole density, n^. In this case, the 
cross section of the Fermi surface A(k z ) has an analyti- 
cal form. The number of electrons per cm 3 is given by 
n e = 4 x 10 24 x fu/v uc with f u — v e /vt z where v e is the 
electron pocket volume and vi, z the BZ volume. v uc is 
the unit cell volume in A~ 3 . Similarly, by replacing v e 
with Vh (the volume of the whole pocket) one can obtain 
nh the number of holes per cm 3 . The critical quanti- 
ties are v e and Vh and they are obtained by integrat- 
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FIG. 7: TTS-GW QP band structure around (a,b) K point, (c,d) H point. The points 1,2 and 3 denote the dispersion we use to 
determine electron and hole masses (see Fig. [9j. The points 1,2 and 3 span a parabola. In (a,b) the parabola is along the KM 
direction with heaviest electron masses and in (c,d) the hole masses are isotropic around H. Note that the parabolas along 
1 — 3 are used for evaluation of the effective electron and hole masses. 



ing the cross section of the Fermi surface, A(k z ) along 
k z . The analytical expression and their dependence on 
the TB parameters is given in [21] . This yields plasmon 
frequencies of huj a = h^J n e e 2 / (m*eo£ a ) = 113 meV for 
plasmon oscillation parallel to the graphcnc layers and 
huj c = hyj n e e 2 / (m* e e e c ) = 19 meV for plasmon oscil- 
lation perpendicular to the layers. Here e a = 5.4 [32] 
and e c = 1.25 [33j are adopted for the dielectric con- 
stants parallel and perpendicular to the graphene layers, 
respectively. We have made two simplifications: first we 
do not consider a finite value of temperature (w and u> c 
are the plasmon frequencies at K) and second we used 
an effective mass averaged over the whole k z range of the 
pockets as shown in Fig(8]and in Fig|9] 

Next we discuss the doping dependence of the elec- 
tronic properties in the so-called dilute limit, which refers 
to a very low ratio of dopant/carbon atoms. Here we use 
a method described previously |21j employing the TB- 
GW parameters from Table [TT] This allows us to calcu- 
late the doping dependence of n e and nh- The analytical 
formula for the Ep dependent cross section of the Fermi 



surface, A(k z ) (given in j^TJ) is integrated for different 
values of Ep. By integrating dA(k z )/dEp along k z we 
obtain the carrier density per eV. The ratio of dopant 
to carbon is given by r — n x v uc /Af where n = n e for 
electron doping and n — for hole doping. Here / is 
the charge transfer value per dopant atom to Carbon. 
Although there are some discussions about the value of 
/, it is was recently found for potassium doping that 
/ = 1 [34J. In Fig. 10 'a) we show the doping dependence 
of the carrier densities. It is clear that at Ep = ±25 meV, 
we have a discontinuity in the carrier density and this is 
associated to the Ep at which the electron or hole pocket 
is completely filled. Since the density of states decreases 
suddenly after the pockets are filled, the kink in the den- 
sity of states appears. This also marks the transition 
from a two-carrier regime to a single-carrier regime. In 
Fig. 10 'b) we show n e and nt- It is clear that at Ep = 
the number of holes equals the number of electrons. At 
Ep = 25 meV, we have no more holes and thus a transi- 
tion from a semi metal to a metal occurs. In such a metal, 
the carriers are electrons, hence an N-type metal. Simi- 
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FIG. 8: Evaluation of the in-plane (a) electron and (b) hole masses around K and H, respectively. The dispersions along the 
parabolas depicted by points 1 — 2 — 3 in Fig. [7] are taken for evaluation of the masses. The effective masses perpendicular to 
the layers in z direction are shown in (c) for electrons and in (d) for holes. □ correspond to TB-GW values and the red lines 
are parabolic fits. 




tions are important for ambipolar transport in graphite 
and graphene: they determine the region for the gate 
voltage in which ambipolar transport is possible. Finally 
in Fig. 10 'c) we plot the ratio of dopant to carbon atoms 



as a function of Ep . 



FIG. 9: The k z dependence of the in-plane electron mass 
(m*) and hole mass (m£) calculated by TB-GW. The masses 
are evaluated along the parabolas as shown in Fig. [8] 



larly, at Ep = —25 meV we have no more electrons and 
the a semi metal to metal transition occurs in the other 
direction to Ep. For this metal, the carriers are holes, 
hence a P-type metal. These semi metal to metal transi- 



It is certainly interesting to monitor the doping in- 
duced changes in n e and also in the shape of the Fermi 
surface. In FigjTTJwe show the Fermi surfaces for electron 
doping and hole doping for Ep = to ±35 meV in steps 
of 5 meV. It is clear that at ±25 meV a single carrier 
regime dominates as indicated by the two different colors 
(red for electrons and green for holes). This is consis- 
tent with Fig. 10 'b) where the integrated electron (hole) 
densities disappear at Ep = 25 meV (Ep = —25 meV). 
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FIG. 10: Doping dependence of (a) the electron and hole carrier density (b) the number of electrons and holes and (c) the 
stochiometric dopant to carbon ratio. The red(green) lines represent electron(hole) carriers. 
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FIG. 11: The Fermi surface of doped graphite in the dilute limit for electron doping and hole doping. The red (green) surfaces 
are the electron(hole) pockets of the Fermi surface. It is clear that at ±25 meV we have a semi metal to metal transition. 



FEW LAYER GRAPHENE 



The 3NN TB-GW set of parameters fits the whole k z 
range of the 3D graphite BZ. Thus the set can be trans- 
ferred for the calculation of QP dispersions of stacked sp 2 
FLGs with N layers (N = 1,2, ... ). The transferability 
of TB parameters is a result of the fact that the lattice pa- 
rameters of FLGs and graphite are almost identical [15]. 
We can use the matrix elements shown in Fig. [lja) also 
for FLGa; for N = 1 only 7q — 7q, sq — s 3 and E are 
needed and this results in the graphene monolayer case. 
For N — 2 the parameters 72 and 75 are not needed 
since they describe next nearest neighbour interactions 
which do not exist in a bilayer. We use the set of TB 
parameters given in Table |T] and the Hamiltonians given 
in section 2 and the appendix. In Fig. [12] we show the 
bilayer (N — 2), the trilayer (N — 3) and the quad- 
layer (N = 4) calculated with TB-GW. The Fig. [L^a) 
and (b) shows the electron dispersion of the bilayer. The 
separation between CB (VB) to the Fermi level is pro- 
portional to 71 [37j and thus it is clear that the TB-GW 
also gives an about 20% larger separation between the 
VB than LDA. Furthermore the slope of all bands be- 
comes steeper for TB-GW since 70 — 7q increase with 
respect to LDA calculations. This is also responsible for 
the increase in vf of FLGs (similar to the graphite case, 
when going from LDA to GW). The same argument is 
the case for the tri- and the quadlayer. It is interesting 
that the QP dispersion measured by ARPES [6] are in 
better agreement with the TB-GW rather than the LDA 



calculations performed. 

The low energy dispersion relation of FLGs are par- 
ticularly important for describing transport properties. 
From the calculations, we find that all FLG has a finite 
density of states at Ep. The trilayer has a small over- 
lap at K point between valence and conduction band (i.e. 
semi-metallic) . This property might be useful for devices 
with ambipolar transport properties. Our results are in 
qualitative agreement with the LDA calculations [48 . 

It can be seen that a linear (Dirac-like) band appears 
for the trilayer (and also for all other odd- numbered mul- 
tilayers). This observation is in agreement to previous 
calculations and is relevant to a increase in the orbital 
contribution to diamagnetism |49j . 

Since TB allows for rapid calculation of the QP bands, 
the transition from FLG to bulk graphite can be ana- 
lyzed. Even in the case of N = 30, the solution of Eq. [4] 
for the 60 x 60 Hamiltonian and overlap matrices takes 
only ~ 10 sec on a Pentium III workstation per k point. 
In Fig. [13] we show the eigenvalue spectrum for FLG with 
N = 1, . . . , 30. As we increase N and hence the number 
of 7T bands, the bandwidth also increases and approaches 
that of bulk graphite. It can be seen that for N > 15, the 
total bandwidth is that of bulk graphite. Interestingly, 
the energies of the tt bands group together and form fam- 
ilies of the highest, second highest etc. energy eigenvalue 
at K. With increasing number of layers, a given family 
approaches the limit for bulk graphite. Such a family pat- 
tern is a direct consequence of the AB stacking sequence 
in FLG and it might be accessible to optical spectroscopy 
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Parameter 


Symbol 


TB-GW 


Experimental value(s) 


Fermi velocity at H [10 6 ms _1 ] 


VF 


1.01 


0.91 :(5 . 1.06 10 . 1.07 [36, 1.02 [37] 


Splitting of 7r bands at K [eV] 


5 


0.704 


0.71 [TO] 


Bottom of 7r band at A point [eV] 


E(A) 


7.6 


8 [38], 8 [39] 


in-plane electron mass [mo] 


m* 


0.1 (fe z averaged) 


0.084 [40], 0.42 [35],0.028 [36] 


in-plane hole mass [mo] a 


ml 


0.06 (k z averaged) 


0.069 05], 0.03 gU 6 , 0.028 [15] 


out-of-plane electron mass [mo] 


m* e 


16 




out-of-plane hole mass [mo] 


™*zh 


-16 




number of electrons at Ef=0 [10 18 cm -3 ] 


n e 


5.0 


8.0 05], 3.1 g2] 


number of holes at E F =0 [10 18 cm" 3 ] 


rih 


5.0 


3.1 [35], 2.7 [42], 9.2 [38] 


Gap at H point [meV] 


A 


5 


5-8 [MEZ] 


in-plane plasmon frequency [meV] 


hiO a 


113 


128 [43] 


out-of-plane plasmon frequency [meV] 


hjjj c 


19 


45-50 [31 HUES] 


A\ optical transition energy [eV] 


Eai 


0.669 


0.722 gS]" 


A-z optical transition energy [eV] 


Ea2 


0.847 


0.926 gg 6 



"To compare different notations, we denote here the absoiute vatue 
of ml . 

Mass at the H point was measured. Note that this is in excellent 
agreement to our calculated k z dependence of mt (see Fig. 9). 

c The difference between the experimental and the T3-GW value 
might be a result of excitonic effects. 

TABLE IV: Properties of the electronic band structure of graphite calculated from TB-GW and compared to experiment. 



similar to the fine structure around K point that has been 
observed in bulk graphite 05] . 



8. DISCUSSION 

We first discuss the QP band structure and relation to 
recent ARPES experiments. From several experimental 
works it is clear that the LDA bands need to be scaled in 
order to fit the experiments [HI EH EE]- Our new set of TB 
parameters quantitatively describes the QP dispersions 
of graphite and FLG. The scaling is mainly reflected in an 
increase of 7q — 7q and 71 , the in-plane and out-of-plane 
coupling, respectively. It thus can be used to analyze 
ARPES of both pristine and doped (dilute limit) graphite 
and FLG. Most importantly, the correlation effects in- 
crease vp, the Fermi velocity when going from LDA to 
GW. For the TB-GW calculation, v F = 1.01 x 10 6 ms" 1 
which is in perfect agreement to the values from ARPES 
that is equal to vp — 1.06 x 10 6 ms _1 [10] . The question 
why LDA works for some metals but fails to give the 
correct vp and energy band dispersion in semi metal- 
lic graphite arises. In graphite, the contribution of the 
electron-electron interaction to the self-energy is unusu- 
ally large. The reason for this is the small number of 
free carriers to screen efficiently the Coulomb interac- 
tion. In most other metals, the density of states at Ep 
has a much larger (~ 1000 times) value than in graphite 
and the screening lengths are shorter. Hence the LDA 
is a good description for such a material but it fails in 
the case of graphite. Another parameter that illustrates 
the quality of the TB-GW to reproduce the experimental 



QP dispersion is 5, the band splitting at K point. For 
TB-GW, we obtain 6 = 0.704 eV and the ARPES gives 
5 = 0.71 eV [H]. 

The TB parameters of the SWMC Hamiltonian have 
been fitted in order to reproduce double-resonance Ra- 
man spectra [14]. As a result they obtained that 73 
(one of the parameters that couples the neighbouring 
graphene planes) has a value of 0.1 eV while the for 
graphite that we obtain in this paper is 0.274 eV (see Ta- 
ble II) and the value that was fit to transport experiments 
is 0.315 eV [THUSO] which is also in perfect agreement to 
ARPES experiments of graphite single crystals [50] . We 
now discuss a possible reason for this discrepancy. The 
fitting procedure in Ref [14] depends on the choice of the 
phonon dispersion relation of graphite. It is important 
to note that recently a Kohn anomaly has been directly 
observed by inelastic x-ray scattering experiments using 
synchrotron radiation |51j which has a steeper slope of 
the TO phonon branch at K point in contrast to pre- 
vious measurements [52]. Thus the assumption of the 
correct phonon dispersion relation is crucial for obtain- 
ing the correct band structure parameters. 

Next we discuss the present QP electronic energy band 
structure in relation to optical spectroscopies such as op- 
tical absorption spectroscopy (OAS) and resonance Ra- 
man spectroscopy. The optical spectroscopies probe the 
joint density of states (JDOS) weighted with the dipolc 
matrix elements. Peaks in the OAS are redshifted when 
compared to the JDOS of the QP dispersion if excitons 
are created. By comparing the QP dispersion with OAS 
experiments [45] [53] [54] we now estimate a value for 
exciton binding energies in graphite. In general many 
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resonant states with k contribute to OAS but a fine 
structure in OAS of bulk graphite measured in reflection 
geometry was observed by Misu et al. [45 which was 
assigned to two specific transitions around K: the A\ 
transition between the lower VB and states just above 
Ep and the A2 transition between the upper VB and 
the upper CB. The experimental energies they found 
were E^f = 0.669 eV and E e jg = 0.847 eV. The TB- 
GW dispersion yields energies E'ffi — 0.722 eV and 
E%Y = 0.926 eV [see Figj^b)]. Assuming one exci- 



ton is created for the A1(A2) transition, this yields ex- 
citon binding energies of — E^f ~ 50 meV and 

E A2 ~ E A2 ~ 80 mcV for a K P° int exciton in bulk 
graphite. Such a value for the exciton binding energies 

most probably increases when going from bulk graphite 

to FLGs due to confinement of the exciton wave function 

in z direction. 

Concerning the value of A, the gap at H, several ex- 
perimental values exist and magneto reflectance exper- 
iments suggest A=5 meV. This is in disagreement to 
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FIG. 13: Evolution of the eigenvalue spectrum for FLG from 
N = 1 ... 30 graphene layers calculated by TB-GW at the 
K point. The dashed lines labelled by iri and 714 denote the 
lower and upper limits of the total bandwidth at K point in 
3D graphite. A pattern that connects the first, second etc. 
energies of FLG is emerging (see text). 



the calculated values obtained for pristine graphite [TP] , 
However, when the doping level is slightly increased, A 
becomes smaller and we thus fixed a doping level in the 
ab-initio calculation that reproduces the experimental 
A [3T] . While some of the variations may be explained 
by the sample crystallinity in z direction, it is also con- 
ceivable that small impurities are responsible for the dis- 
crepancy. 

Next we discuss a possibility to measure the free car- 
rier plasmon frequencies of the of pristine and alkali- 
metal doped graphite by high resolution energy electron 
loss spectroscopy (HREELS). The electron concentration 
inside the pockets very sensitively affects the plasmon 
frequencies. In principle the charge carrier plasmons 
should also appear as a dip in optical reflectivity mea- 
surements but due to the small relative change in in- 
tensity they have not been observed so far. Due to the 
small size of the pockets, the number of charge carriers 
and hence the conductivity and plasmon frequencies are 
extremely sensitive to temperature and doping. Exper- 
imentally observed plasmon frequencies for oscillations 
parallel to the graphene layers are Sw a =128 meV [33J and 
for oscillations perpendicular to the layers are hw c =Ab — 
50 meV (33J HIJ [44] . These values for hio a and Hu> c agree 
reasonably well with our derived value considering that 
we make the crude estimation of huj a and huj c at K while 
experiments where carried out at room temperature. We 
also used average effective masses of the electron pocket 
in order to determine the plasmon frequency, while in 
fact there is a k z dependence (see Figj9| . We also note 
that the experimental literature values for e a and e c have 
a rather wide range, e.g. e c = 3.4 |55] from reflectivity 



measurements and e c = 5.4 [32] from EELS. 

The temperature dependence of uj c was measured by 
Jensen et al [33 by HREELS and they observed a strong 
T dependence for lo c which was attributed to changes in 
the occupation and thus number of free carriers with T. 
The observed plasmon energy was rising from 40 meV to 
100 meV in a temperature range of 100 K to 400 K. A 
similar effect might be observed by HREELS of doped 
graphite as a function of doping level. It would be inter- 
esting to study the evolution of the plasmon frequency 
with doping level. With our current understanding of 
the low-energy band structure we predict a semimctal to 
metal transition at a Fermi level shift of Ep ~ 25 meV. 
At this doping level, the hole pocket is completely filled 
with electrons and disappears and the electron pocket 
has roughly doubled in size and u> a and uj c should in- 
crease by a factor of ~ V2- Such a transition should be 
observable by HREELS and might be accompied by inter- 
esting changes in the band structure (electron-plasmon 
coupling) which can also be measured simultaneously by 
ARPES. Many DC transport properties can be under- 
stood with a Drude model for the conductivity, which 
is inversely proportional to the effective carrier mass. 
From the ratio of m* e /m* and TO* h /m^, we hence expect 
that the DC electron(hole) conductivity in z direction is 
~ 200(~ 500) times less than the in-plane conductivity 
. The experimental value is 3 x 10 3 [19] and considering 
the large variation in experimental values even for the in- 
plane conductivity |19M56j . it is in reasonable agreement. 
In Table II VI we summarize the electronic band structure 
properties such as Vp, 5, the effective masses etc. and 
compare them to experimental values. 



9. CONCLUSIONS 

In conclusion, we have fitted two sets of TB parameters 
to first principles calculations for the bare band (LDA) 
and QP (GW) dispersions of graphite. We have observed 
a 20% increase of the nearest neighbour in-plane and 
out-of-plane matrix elements when going from LDA to 
GW. Comparision to ARPES and transport measure- 
ments suggest that LDA is not a good description of the 
electronic band structure of graphite because of the im- 
portance of correlation effects. We have explicitely shown 
that the accuracy of the TB-GW parameters is sufficient 
to reproduce the QP band structure with an accuracy 
of ~ 10 meV for the higher energy points and ~ 1 meV 
for the bands that are relevant for electronic transport. 
Thus the 3NN TB-GW calculated QP band energy dis- 
persions are sufficiently accurate to be compared to a 
wide range of optical and transport experiments. The 
TB-GW parameters from Table I should be used in the 
future for interpretation of experiments that probe the 
band structure in graphite and FLG. For transport ex- 
periments that only involve electronic states close to the 
KH axis (in graphite) or the K point in FLG, the SWMC 
Hamiltonian with parameters from Table II can be used. 
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With the new set of TB parameters we have calculated 
the low energy properties of graphite: the Fermi velocity, 
the Fermi surface, plasmon frequencies and the ratio of 
conductivities parallel and perpendicular to the graphcnc 
layers and effective masses. These values are compared 
to experiments and we have demonstrated an excellent 
agreement with almost all values. 

We have shown that at ±25 meV, a semi metal to metal 
transition exists in graphite. Such a transition should 
be experimentally observable by HREELS or ARPES for 
electron doped graphite single crystals with very low dop- 
ing levels (i.e. the dilute limit). Experimentally, the 
synthesis of such samples is possible by evaporation of 
potassium onto the sample surface and choosing a suffi- 
ciently high equilibration temperature so that no staging 
compound forms. It would thus be interesting to do a 
combined ARPES and HREELS experiment on this kind 
of sample. 

We have also discussed the issue of the existence of 
Dirac Fermions in graphite and FLGs. For pristine 
graphite it is clear that for a finite value of the gap A at 
H point, we always have a non-zero value of the effec- 
tive masses at H. We now discuss some possibilities for 
A = 0. One way would be to use AA stacked graphite. 
However, natural single crystal graphite has mainly AB 
stacking and the preparation of an AA like surface by 
repeated cleaving is a rather difficult task. Recently, the 
possibility of reducing A by increasing the doping level 
has been put forward [3TJ [34] and this seems a more 
promising way for experiments. For the case of a fully 
doped graphite single crystal, the interlayer interaction 
is negligible because the distance between the graphene 
sheets is increased and the stacking sequence is change 
to AA stacking. These two effects would cause the ap- 
pearance of Dirac Fermions in doped graphite. 

With the new set of 3NN TB parameters we have calcu- 
lated the QP dispersion of FLG in the whole 2D BZ. We 
have shown the evolution of the energy eigenvalue spec- 
trum at K point as a function of N, the number of layers. 
An interesting family pattern was observed, where the 
first, second, third etc. highest transitions form a pattern 
that approaches the energy band width of bulk graphite 
with increasing N. For the highest transition, the bulk 
graphite band width is already reached at N = 15. This 
family pattern is a direct result of the AB stacking se- 
quence in FLG and it can possibly be observe experi- 
mentally. It should be pointed out here, that a very sim- 
ilar fine structure around the K point in bulk graphite 
has been observed by optical absorption spectroscopy by 
Misu et al. [IS] although the family pattern from Fig. 13 
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might be observed also by ARPES on high quality FLGs 
synthesized on Ni(lll). 
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Appendix 

We made computer programs to automatically gener- 
ate the 3NN Hamiltonians and overlap matrices for N 
graphene layers. The method is following the descrip- 
tion by Partoens et al. [TB] but with the difference that 
we include 3NN matrix elements that describe the QP 
dispersion in the whole BZ. It is interesting to note the 
similarity between the Hamiltonian of a N layer graphene 
and a stage N graphite intercalation compound |5Sj. The 
general form of the Hamiltonian is given by 



/ Ha 1 a 1 Ha 1 b 1 Hau 



ff(k) 



H 



BiAi 



H 



B 1 B t 



Ha 2 Ai Ha~ 



Si 



Hb!A 2 

Ha 2 a 2 



V-fffiivAi H Bn B! H Bn a 2 
and the overlap matrix is given by 



5(k) 



Sa 1 A 1 SU1B1 SaxA 2 

SbxAi Sb 1 B 1 Sb 1 A 2 

Sa 2 Ai Sa 2 Bi Sa 2 a 2 

\ Sb n A! S Bn b 1 S Bn a^ 



HaxB n \ 

Hb 2 b n 
Ha 2 b n 

Hb n Bn / 



Sa!B n \ 
S B 2 B N 
Sa 2 b n 

S Bn Bn ) 



tion Eq. [4] The 3NN TB Hamiltonians can be used with 
GW parameters from Table |T] In the following we give 
the explicit form for the H(k) and S(k) of a monolayer 
and an AB stacked bi- and and ABA trilayer which was 
used in section 13 



Monolayer graphene 

The monolayer 3NN Hamiltonian and overlap matrices 
are given by 



^o+7o 2 / 2 (k) 7o 1 /i(k)+7o 3 / 3 (k) 
7 1 / 1 (k)*+ 7o 3 / 3 (k)* £ + 7o 7 2 (k)* 




l+/ 2 (k)«g 



/i(k) a J+/ 3 (k) s g 



/i(k)* s J+/ 3 (k)* a g l+/ 2 (k)* S g 



The QP band structure is then given by solving equa- 



Bilayer graphene 

The bilayer 3NN Hamiltonian and overlap matrices are 
given by 



H(k) = 

( £ + A+7 2 / 2 (k) 
7o 1 /i(k)*+7o 3 / 3 (k)* 

7i 

V 74/i (k) 



7o7i(k)+7o 3 / 3 (k) 

i?o+7o 2 / 2 (k)* 
74/i (k) 
73/1 (k)* 



7i 74/1 (k)* \ 

74/i(k)* 73/1 (k) 

E + A+ 7o 2 / 2 (k) 7o 1 / 1 (k)*+ 7o 3 / 3 (k)* 
7 1 / 1 (k)+ 7o 3 / 3 (k) So+7o 2 / 2 (k) / 



(5) 



and 



S(k) 
( 



l+sg/ 2 (k) 4/ 1 (k)+ S 3/ 3 ( k ) 



V 



aS/i(k)*+ S g/ 3 (k)* l+*g/ 2 (k)* 





4/ 1 (k)+ s 3/ 3 (k) i+4f 2 (k) J 



(6) 



1G 



Trilayer graphene 

The trilayer Hamiltonian and overlap matrices read as 
H(k) = 

( E + A+ 7o 2 /2(k) / 1 (k) 7o 1 + 7o 3 / 3 (k) 71 74/i (k)* 75 \ 

7o 1 A(k)*+7 3 / 3 (k)* i?o+ 7o 2 /2(k)* 74/1 (k)* 73/1 (k) 72 

7i 74.fi (k) £o + A+7 2 / 2 (k) 7o 1 /i(k)*+7o 3 / 3 (k)* 71 7 4 /i(k) 

74/i (k) 73/1 (k)* 7o 1 /i(k)+7o 3 /3(k) S +7o 2 / 2 (k)* 74/1 (k) 73/1 (k)* 

75 71 74/i (k)* £ + A+ 7o 2 / 2 (k) 7o 1 /i(k)+7o 3 /3(k) 

V 72 74/i (k)* 73/i(k) 7o 1 /i(k)*+ 7o 3 / 3 (k)* £ 0+7 2 / 2 (k)* J 



and 



S(k) = 

/ l+*g/ 2 (k) /i(k)«5+ag/ 3 (k) 

*J/i(k)*+ag/ 3 (k)* l+«g/ 2 (k)* 

l+*§/ 2 (k) 4.f 1 (k)*+ s 3 / 3 (k)* 

*J/i(k)+ag/ 3 (k) l+^/ 2 (k) 

£ + A+s 2 / 2 (k)* *J/i(k)+ag/ 3 (k) 

V *S/i(k)*+*§/ 3 (k)* l+«§/ 2 (k)* 







The sum of the phase factors for first-, second- and 
third nearest neighbours are given by /i(k), / 2 (k) and 
/ 3 (k), respectively. Note that /i(k) and / 3 (k) couple 
A and B atoms and / 2 (k) describes the AA and BB 
interactions. 



/ 3 (k) = ^exp(ikr|) 



(11) 



/.i.(k) = oxp ( i ]+2rxp 



ik x a a \ ( y/3k y a 
cos -r — 



/ 2 (k) = ]Texp(zkr 2 ) 



(9) 



(10) 



Here r 2 and r| are the vectors that connect the A\ 
atom (see Fig. 1) with the second nearest B atoms (three 
atoms) and the third nearest A atoms (six atoms), re- 
spectively. 



[1] Th. Seyller, K.V. Emtsev, K. Gao, F. Speck, L. Ley, 
A. Tadich, L. Broekman, J.D. Riley, R.C.G. Leckey, 
O. Rader, A. Varykhalov, and A.M. Shikin, Surf. Sci. 
600, 3906 (2006). 

[2] A. Griineis and D. Vyalikh, Phys.Rev. B 77, 193401 
(2008). 

[3] K.S. Novoselov, D. Jiang, F. Schedin, T.J. Booth, V.V. 
Khotkevich, S.V. Morozov, and A.K. Geim, PNAS 102, 
10451 (2005). 

[4] A. Geim and K. Novoselov, Nature Mat. 6, 183 (2007). 

[5] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Na- 
ture Phys. 2, 620 (2006). 

[6] T. Ohta, A. Bostwick, T. Seyller, K. Horn, and E. Roten- 
berg, Science 313, 951 (2006). 

[7] C. Berger, Z. Song, T. Li, X. Li, A.Y. Ogbazghi, R Feng, 



Z. Dai, A.N. Marchenkov, E.H. Conrad, P.N. First, and 
W.A. de Heer, J. Phys. Chem. 108, 19912 (2004). 
[8] V. M. Karpan, G. Giovannetti, P. A. Khomyakov, M. Ta- 
lanana, A. A. Starikov, M. Zwierzycki, J. van den Brink, 

G. Brocks, and P. J. Kelly, Physical Review Letters 
99(17), 176602 (2007). 

[9] A. Bostwick, T. Ohta, T. Seyller, K. Horn, and E. Roten- 
berg, Nature Phys. 3, 36 (2007). 
[10] A. Griineis, C. Attaccalite, T. Pichler, V. Zabolotnyy, 

H. Shiozawa, S.L. Molodtsov, D. Inosov, A. Koitzsch, 
M. Knupfer, J. Schiessling, R. Follath, R. Weber, 
P. Rudolf, L. Wirtz, and A. Rubio, Phys. Rev. Lett. 100, 
037601 (2008). 

[11] S. Y. Zhou, G. H. Gweon, C. D. Spataru, J. Graf, D. H. 
Lee, S. G. Louie, and A. Lanzara, Phys. Rev. B 71, 



17 



161403 (2005). 

[12] K. Sugawara, T. Sato, S. Souma, T. Takahashi, and 

H. Suematsu, Phys. Rev. Lett. 98, 036801 (2007). 
[13] S. Reich, C. Thomsen, and P. Ordejon, (2002). Phys. 

Rev. B 65, 153407 (2002). 
[14] L.M. Malard, J.Nilsson, D.C.Elias, J.C.Brant, F.Plentz, 

E.S.Alves, A.H.Castro Neto, and M.A.Pimenta, Phys. 

Rev. B 76, 201401 (2007). 
[15] A. Bostwick, T. Ohta, J.L. McChesney, T. Seyller, 

K. Horn, and E. Rotenberg, Solid State Communications 

143(1-2), 63-71 (2007). 
[16] B. Partoens and F.M.Peeters, Physical Review B 74, 

75404 (2006). 

[17] K. Sasaki, J. Jiang, R. Saito, S. Onari, and Y. Tanaka, 

J. Phys. Soc. Jpn. 76, 033702 (2007). 
[18] R.C. Tatar and S. Rabi, Physical Review B 25, 4126 

(1982). 

[19] M. S Dresselhaus and G. Dresselhaus, Advances in Phys. 
30, 139 (1981). 

[20] M.S. Dresselhaus and J. Mavroides, IBM Journal of Re- 
search and development page 262 (1964). 
[21] M. S. Dresselhaus, G. Dresselhaus, and J. E. Fischer, 

Phys. Rev. B 15, 3180-3192 (1977). 
[22] J.C. Charlier, X. Gonze, and J. P. Michenaud, Phys. Rev. 

B 43, 4579 (1991). 
[23] S. Y. Zhou, G. H. Gweon, J. Graf, A. V. Fedorov, C. D. 

Spataru, R. D. Diehl, Y. Kopelevich, D. H. Lee, S. G. 

Louie, and A. Lanzara, Nature Phys. 69, 245419 (2006). 
[24] W. W. Toy, M. S. Dresselhaus, and G. Dresselhaus, Phys. 

Rev. B 15, 4077-4090 (1977). 
[25] G. P. Mikitik and Yu. V. Sharlai, Phys.Rev. B73, 235112 

(2006). 

[26] X. Gonze, J.M. Beuken, R. Caracas, F. Detraux, 
M. Fuchs, G.M. Rignanese, L. Sindic, M. Verstraete, 

G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami, Ph. 
Ghosez, JY. Raty, and D.C. Allan, Comput. Mater. Sci. 
page 478 (2002). 

[27] M .S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 
(1986). 

[28] L. Hedin, Phys. Rev. 139, A796-A823 (1965). 

[29] S. G. Louie, Topics in Computational Materials Sci- 
ence, edited by C. Y. Fong (World Scientific, Singapore) 
page 96 (1997). 

[30] A. Marini et al., the Yambo project, |http: / /www.yambo- 
code. org/| 

[31] C. Attaccalite, A. Griineis, T. Pichler, and A. Rubio, in 

preparation (2008). 
[32] H. Venghaus, Phys.Stat.Sol.B 71, 609 (1975). 
[33] E.T. Jensen, R.E. Palmer, and W. Allison, Phys. Rev. 

Lett. 66, 492 (1991). 
[34] A. Griineis, C. Attaccalite, A. Rubio, D. Vyalikh, 

S.L. Molodtsov, J. Fink, R. Follath, W. Eberhardt, 

B. Biichner, and T. Pichler, Nature Mat. in review (). 
[35] S. Y. Zhou, G. H. Gweon, J. Graf, A. V. Fedorov, C. D. 

Spataru, R. D. Diehl, Y. Kopelevich, D. H. Lee, S. G. 

Louie, and A. Lanzara, Nature Phys. 69, 245419 (2006). 
[36] G. Li and E. Andrei, Nature Phys. 3, 623 (2007). 
[37] M. Orlita, C. Faugeras, G. Martinez, D.K. Maude, M.L. 

Sadowski, and M. Potemski, Phys. Rev. Lett. 100, 

136403 (2008). 

[38] K. Sugawara, T. Sato, S. Souma, T. Takahashi, and 

H. Suematsu, Physical Review B 73, 045124 (2006). 
[39] A.R. Law, M.T Johnson, and H.P. Hughes, Physical Re- 
view B 34, 4289 (1986). 



[40] E. Mendez, T.C. Chieu, and N. Kambe M.S. Dresselhaus, 

Solid State Commun. 33, 837 (1980). 
[41] J. K. Gait, W. A. Yager, and H. W. Dail, Phys. Rev. 

103(5), 1586-1587 (Sep 1956). 
[42] D.E. Soule, Phys. Rev. 112, 698 (1958). 
[43] J. Geiger, H. Katterwe, and B. Schroder, Z.Phys.B 241, 

45 (1971). 

[44] P. Laitenberer and R. E. Palmer, Phys. Rev. Lett. 76, 
1952 (1996). 

[45] E. Mendez, A. Misu, and M. S. Dresselhaus, Phys. Rev. 

B 21, 827 (1980). 
[46] S.B. Trickey, F.Miiller-Plathe, G.H.F. Diercksen, and 

J.C. Boettger, Physical Review B 45, 4460 (1992). 
[47] Edward McCann and Vladimir I. Falko, Phys. Rev. Lett. 

96, 086805 (2006). 
[48] Sylvain Latil and Luc Henrard, Phys. Rev. Lett. 97, 

036803 (2006). 

[49] Mikito Koshino and Tsuneya Ando, Physica E: Low- 
dimensional Systems and Nanostructures 40(5), 1014 
(2008). 

[50] A. Griineis, T. Pichler, H. Shiozawa, C. Attaccalite, 
L. Wirtz, S.L. Molodtsov, R. Follath, R. Weber, and 
A. Rubio, physica status solidi (b) 244, 41294133 (2007). 

[51] A. Griineis et al., submitted (2008). 

[52] J. Maultzsch, S. Reich, C. Thomsen, H. Requardt, and 
P. Ordeon, Phys. Rev. Lett. 92, 75501 (2004). 

[53] A. Misu, E. Mendez, M. S. Dresselhaus, and K. Nakao, 
J. Phys. Soc. Jpn. 47, 199-216 (1979). 

[54] J. M. Zhang and Eklund, J.Mat.Res. 2, 858 (1987). 

[55] M. Zanini and J. Fischer, Mat.Sci.Eng. 31, 169 (1977). 

[56] D.T. Morelli and C. Uher, Phys. Rev. B 30, 1080 (1984). 

[57] A. Kokalj, Comp. Mat. Sci. 28, 155 (2003). 

[58] R. Saito and H. Kamimura, Phys. Rev. B 33(10), 7218 
(1986). 



